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ABSTRACT 



This paper presents a real time seismic event detection 
and source location (RSEDSL) algorithm using single station 
(three component) data for seismic analysis, RSEDSL detects 
and locates an event through parameters derived from the data 
covariance matrix in a manner analogous to beam steering 
using data from an array of single component stations. The 
algorithm has been implemented on a TMS32020 real time 
digital signal processing (DSP) system and tested on several 
data sets. The algorithm has also been implemented on a Sun 
2/170 mainframe. Comparison of the results obtained fiom the 
TMS with those obtained from the Sun indicate that seismic 
events can be accurately detected and located in real time. 
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1. Introduction 

Seismic analysis using single station (three component) 
seismic data has important appUcations as a means of verifying 
a nuclear test ban treaty (1J. The objective of this paper is to 
show that seismic events can be accurately detected and located 
in real time. To accomplish this objective RSEDSL has been 
implemented on a TMS32020 [2.3] DSP system and tested on 
several seismic events. The results were then compared to 
those obtained from an equivalent implementation on a Sun 
2/170 mainframe. The results from the the real time imple- 
mentation on the TMS32O20 were comparable to those 
obtained from the Sun. 

RSEDSL is designed to detect and locate regional seismic 
events using short period seismic data (sampled at 40 Hz) as 
input. A regional seismic event is de5ned to be one in which 
the angle subtended at the center of the earth by an arc con- 
necting the receiver and source (epicenter) is less than 20°. 
This angle (4) corresponds to a surface distance of 2224 km 
[41. The algorithm has been extensively tested using data 



obtained from a regional seismic test network (RSTN) origi- 
nally designed and built by Sandia National Laboratories 
(Albuquerque, NM 87185). Each station in this network 
senses earth motion along three orthogonal axes North (N). 
East (E), and vertical (Z). This data was originally sampled at 
40 Hz and later lowpass altered at 3.5 Hz and resampled at 8 
samples per second (sps). As currently implemented, RSEDSL 
uses only the N and E data channels for event detection. - Once 
an event is detected, the instantaneous source azimuth ($<m)) is 
determined by estimating the polarization direction of the initial 
compression^ phase (/>„). The source azimuth is shown as i 
in Figure I. The Z-channel is used to resolve an inherent 18CT 
ambiguity in fan). The distance to the source is then deter- 
mined from the time difference between the P m phase and a 
later phase (L,). RSEDSL is also capable of detecting inter- 
mediate (S„) phases. 

2. The algorithm 

RSEDSL begins by estimating the input data covariance 
matrix <?(«) of the N and E input channels. C(m) is defined 
as 



C(m)= f*A»> <M«>1 



-On) 6}{m) 

where m is the time ind ex , d 2 (m) u Ae N cfumnel (n(m)) VJU . 
ance. *,V> is the E channel (e<m)} variance, and is the 

N, E channel data crossconelaoon. Note that each of the data 
channels is assumed to be zero mean. The following recursive 
relations are used to compute each element of <?(m) [1] 

«2(m) - - n \m)] + n\ m ) (2a ) 

*frO-PUftm-l)-«V)] + eVi) (2b) 
d ~ (m) " Pt*i(«-D - n(«V(m)] + n(m) e (m) (2c) 
where p is a smoothing parameter with a value between 0 and 
I. For the implementation discussed in this paper 0 has a 
value of 0.98. 
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Event detection is accomplished through the use of a 
detection parameter derived from the maximum eigenvalue 
of C(m). is obtained by solving the charac- 

teristic equation of C(m) as follows 

i Jf {rrT} _ gjfoO + *ftw) 4- VW(m) - ^(m))^ + 4d^Q 

A smoothed delayed value of ^(m) is then computed recur- 
sively as 

**("0 = P[lj(m-1) - ^(m-£>)] + ^(m-D) (4) 
In equation (4) p is a smoothing constant equal to 0.98 and D 
is a delay parameter equal to 200 samples. The detection 
parameter is defined as 

L(rn)=X M ( m )-k k i s (m)>0 ? (5) 
RSEDSL declares an event when k(m) exceeds zero. Once an 
event has been declared, a new event flag is set and detection 
is disabled until returns to its pre-event noise level. The 

pre-event noise level is determined by the value of i,<«) at the 
time the event is declared. 

The onset time of an event is determined by backing up a 
maximum of 10 seconds prior to the detection of the event and 
computing a second detection parameter. This detection 
parameter (&(«)) is computed as in equation (5) above except 
that the scalar multiplier is replaced by a smaller value 
In the current implementation k h is 3.5 and k t is 13. 

After RSEDSL has detected the onset of an event, it then 
obtains an estimate of the direction of travel of the seismic sig- 
nal in the NE-plane as measured from the N axis. The tech- 
nique that is used is analogous to beamforming with a seismic 
array. The instantaneous azimuth is computed from the 

eigenvector corresponding to and is given by 



i(m ) » tan 



(6) 



An estimate of fa) is obtained by averaging equation (6) over 
a period of 2.0 sec. starting 1.0 sec. after the onset time of the 
event [1J. 

As previously mentioned the distance to the source of an 

event is determined from the time difference between the P 

and L t phases. RSEDSL is capable of detecting and 

differentiating these two phases. Standard tables are then used 

to obtain an estimate of the sourcc-to-receiver distance [4] 

™ S ^ has tecn left » *e analyst in the current implemen- 
cation. 



The algorithm has been tuned to detect events with body 
wave magnitudes (m 4 ) i„ die range 3 . 0 w 6 Q Tests using a 
Fortran implementation of the algorithm on a Sun 2/170 com- 
puter have shown that events having a source-to-rccciver dis- 
tance of 1300 km and m> in the range 4.0 to 6.0 could be con- 
sistemly detected and located to within a radius of approxi- 
mately 70 km of the source [1]. The results obtained from the 
implementation on the Sun 2/170 were used as benchmarks 
with which to measure the performance of the real time i mplc . 
mentation on the TMS32020. 

3. Implementation 

RSEDSL was implemented on a 20 MHz TMS32020 

based DSP system supplied by Pacific Mictocircuits Ltd. (6J. 
The system was designed to fit in an expansion slot of an IBM 
PC-XT. Data was transferred between the PC and Sun 2/170 
via an Ethernet link. In addition to the TMS32020's 544 
worts of on-chip RAM. the board is equipped with' I6K of 
high speed external RAM. The board also has 16 bit A/D and 
D/A ports with analog to digital conversion , imc of 17 
Note that for the implementation discussed in this paper the 
input data was 'sampled' from the on-board external RAM and 
not the A/D. 

The most serious problems encountered during the imple- 
mentation of the algorithm were tmncation and roundoff errors 
resulting from the conversion of 32 bit products to 16 bits. 
Initially, these errors were large enough to cause RSEDSL to 
yield unsatisfactory results, especially in events with a weak P. 
phase (7J. Reducing the efTects of these errors required the u« 
of 32 bit arithmetic wherever possible and altering the form of 
the equations used to estimate the covariance matrix. These 
equations (2a-c) were written in the following form )n the ori- 
ginal algorithm 

y(m) = 0y<«-l) + (l-p)xV) (7) 
Implementing this equation in Fortran on the Sun 2/170 
presented no problems since the floating point data format pro- 
vided sufficient dynamic range. However, this is not a very 
efficient form when implemented on the fixed point TMS32020 
and results in serious underflow problems for small input 
values. r 

Whenever the TMS32020 multiplies two 16 bit numbers 
the resulting 32 bit product is saved in a product register 
whose contents can be directly added to the contents of the 32 
bit accumulator. This property was used to reduce underflow 
errors in estimating the covariance matrix elements by implc- 
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mcnting equation (7) in the following form 

y M m PCy(m-l) - x 2 (m)l + x\m) 



(8) 



One can see in this equation that if the output is saved in 32 
bit format, then the term in brackens can be evaluated as a 32 
bit number before being converted to 16 bits for multiplication 
with p. 

Another property of the TMS32020 that helped facilitate 
the use of 32 bit arithmetic in RSEDSL was the normalization 
instruction. This instruction was used to convert 32 bit 
numbers to 16 bit numbers for multiplication by shirting left a 
number of times until a one was in the highest bit position. 
The number of left shifts was saved in an auxiliary register to 
be used for denormalization after multiplication. 

A special approximation subroutine was used to imple- 
ment the square root term in equation (3). This term was of 
the general form 

R = JFTqI (9) 

Many methods exist for obtaining a close approximation of R 
if / and Q are known (8 J. For the implementation discussed in 
this paper / and Q were determined via the following relations 

/ * l*foi)-*ftn)l (10a) 

Q * I2$V (m) I (10b) 

Equation (9) was implemented by using the following formulas 
[8J. 



R |/|>|<2, 

R=\Ql+JLL |<2 1 z> |/| 



(Ua) 
(lib) 

Equations (lla,b) were very simple to implement on the 
TMS32020. Since division by two only requires a right shift 
of one bit, R can be determined by a shift and an addition. 
This approach proved to be better than other schemes tried [7J. 

Computation of #«) in equation (6) also required the use 
of an approximation subroutine to implement the arctangent 
function. This subroutine used one of two approximation tech- 
niques depending upon the magnitude of the argument (*): 

0.0 Sx < 03 : arctan(x) = x 

03 5 x < 1.0 : use look-up table. 250 entries 
resolution = .0028 rad. 
Prior to calling the subroutine the magnitudes of the numerator 
and denominator in equation (6) were compared. If the 
numerator was smaller, then x was less than one and an 
approximation was determined by the appropriate technique as 



indicated above. If the numerator was larger than the denomi- 
nator, then the two were swapped prior to dividing and a flag 
was set An approximation was then made by one of the 
above techniques and the following relation was used to 
account for the swap 

<*=9O*-0 (12) 
Figure 2 shows a plot of the absolute error between the output 
of this subroutine and that of a standard C library arctan func- 
tion evaluated on the Sun 2/170. Note that this error is less 
than # (.00873 rad.). This approximation scheme was found 
to be more accurate and simpler to implement than other 
schemes tried (7). 

The entire program (RSEDSL) was stored in the high 
speed external memory while the arctan look-up table, derived 
parameters, and all temporary variables were stored in the on- 
chip RAM. 

4. Results 

RSEDSL took 220 us to process each input sample vector 
<n(m). e(m). z(ra)) if no event was detected. This figure 
increased to 340 us if an event was detected. This time 
includes determination of onset time and azimuth averaging 
and is independent of which phase was detected (/>„, s m . or 
L g ). These figures indicate that RSEDSL is capable of pro- 
cessing data from a network of several stations even at the ori- 
ginal sampling rate of 40 sps. 

figures 3, 4. and 5 show plots of some of the outputs 
^obtained on the Sun 2/170 and TMS32020. The input used to 
obtain these results corresponds to a nuclear test at the Nevada 
test site observed at an RSTN station in Albuquerque, NM. 
The source-to-receiver distance in this case is approximately 
.800 km. The event had a body wave magnitude of 4.2 [4]. 
Figure 3 shows the input 2-channcl, maximum eigenvalue, 
smoothed delayed maximum eigenvalue, and azimuth estimate 
from the TMS32020 for the entire event. Figure 4 shows the 
same outputs expanded about the P a phase but with the outputs 
from the Sun 2/170 (dotted line) superimposed. Figure 5 also 
shows the same outputs as Figure 3 but expanded about the L t 
phase and with the outputs from the Sun 2/170 (dotted line) 
superimposed. Note that the maximum eigenvalue plots are 
identical, while the azimuth plots differ before the onset of the 
event but are nearly identical once the event starts. Table 1 
shows a comparison of detection times and azimuth estimates 
for each phase of the event as determined on the Sun 2/170 
and TMS32020. Note that detection times on the TMS32020 
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arc within 0.125 seconds and the azimuth estimates arc within 
three degrees for the important P m and L, phases. RSEDSL 
was tested on two additional seismic events (nuclear tests) and 
yielded similar results (7J. 

5. Conclusions 

RSEDSL has been tested on data corresponding to actual 
underground nuclear tests and acheived results comparable to 
those obtained on a Sun 2/170 computer. The execution time 
was short enough to allow RSEDSL to process multiple sta- 
tions in real time. Finally, RSEDSL needs to be tested on a 
larger number of events with varying magnitudes and source- 
to-receiver paths. 
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TMS32020 


SUN 2/170 


P m PHASE 
ONSET TIME 
NE AZM 


85.5 
81.5 
107.51 


85.375 
81.625 
104.49 


NEXT PHASE 
ONSET TIME 
NE AZM 


115375 

105.5 

123.72 


115.5 

105.625 

100.96 


L x PHASE 
ONSET TIME 
NE AZM 


216.0 
214.5 
37.91 


216.125 
214.625 
40.04 



Time in seconds. NE AZM in degrees. 
Table I. Comparison of output. 
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